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^^ We propose a numerical procedure to study closure approximations for FENE dumb- 

^H bells in terms of chosen macroscopic state variables, enabling to test straightforwardly 

which macroscopic state variables should be included to build good closures. The 

method involves the reconstruction of a polymer distribution related to the conditional 

'~^ equilibrium of a microscopic Monte Carlo simulation, conditioned upon the desired 

O . macroscopic state. We describe the procedure in detail, give numerical results for 

several strategies to define the set of macroscopic state variables, and show that the 

resulting closures are related to those obtained by a so-called quasi-equilibrium approx- 



C^ imation [19]. 



> 1 Introduction 

^*>0 The simulation of dilute solutions of polymers in a Newtonian solvent is a challenging mod- 

elling and numerical problem, since deformation of the polymer molecules causes stresses 
(^ that result in macroscopic non-Newtonian rheological behavior. One approach is to couple 

^ the macroscopic fluid flow equations to a microscopic model for the polymers, a so-called 

^1^ micro-macro model [15, 27, 28]. The simplest microscopic models, that we will use in 

this paper, describe the individual polymers as non-interacting dumbbells, consisting of 
H two beads connected by a spring that models intramolecular interaction. The state of the 

cd 

polymer chain is described by the end-to-end vector Xj that connects both beads whose 
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evolution is modelled using a stochastic differential equation (SDE): 
dX.t + u • V^Xi dt 



V^uXt - ^F(Xi) 



dt + J^^dWt, (1.1) 



where u is the velocity field of the solvent, C is a friction coefficient, T is the temperature, 
ks is the Boltzmann constant, and Wt is a standard multidimensional Brownian motion. 
This model takes into account Stokes drag (due to the solvent velocity field), a spring force 
F and Brownian motion (due to collisions with solvent molecules). The left-hand side of 
Equation (1.1) is the convective derivative. Note that the stochastic process Xf implicitly 
depends on the space variable x. 

To specify the microscopic model (1.1) completely, we need to define the spring force. 
This force can be more or less complicated, depending on the effects taken into account. 
The simplest model is the Hookean dumbbell model for which the spring is linear elastic: 

F(X) = iJX, 

with H a spring constant. Another model, which is the focus of this paper and which is 
known to yield better agreement with experiments, is the finitely extensible nonlinear elastic 
(FENE) force [4]: 

^^^^ = i-\\x\\yibkBT/Hy ^^-^^ 

where 6 is a nondimensional parameter related to the maximal polymer length. 

In the macroscopic part of the model, the evolution of the solvent velocity and pressure 
fields u and p is modeled by mass and momentum conservation equations: 

p f — + u • V^^u j = r]s Ar,u - VxP + div^ (Tp), 

^ div^.(u) = 0, 

with p the density and rjg the viscosity. Equation (1.3) contains an additional stress tensor 
Tp due to polymer deformation, which is given via the classical Kramers' expression 

Tp{x,t) = n{-Kt ^F{Xt)) - nkeTld. (1.4) 

Here, n is the polymer concentration and (•) denotes the expectation over configuration 
space, which is approximated in practice by an empirical mean over a very large ensemble 
of realizations of Xt, solutions to (1.1). 

One thus obtains a coupled system (1.1)-(1.3)-(1.4) that we rewrite in a non-dimensional 



form as (see for example [20]): 



Re f ^ + u • V^.u j = (1 - e)A^u - V^p + div^.(Tp), (1.5) 

div(u) = 0, (1.6) 

rp = ^((Xi0F(Xi))-Id), (1.7) 



dX.t + u • Vx'Xi dt 



V^.uXf ^F(Xt) 

2We ^ ' 



dt+—^dWt, (1.8) 

VWe 



where the nondimensional parameters are: 



Re=- , We = ^, e=-^. 1.9 



Here, U is a characteristic velocity, L = y^ksT/H denotes a characteristic length, A = C/'^H 
is a characteristic relaxation time for the polymers and rjp = nksTX is a viscosity associated 
to the polymers. The total viscosity is r] = rjp -\- r]s. The parameters Re and We are the 
Reynolds and Weissenberg number, respectively. The nondimensional Hookean and FENE 
forces write respectively: 

FhookO^) = X, FfeneO^) = 1|„||9 ,, ■ (1-10) 

1 — ||X||^/o 

The microscopic part of the model, i.e. (1.7)-(1.8), can equivalently be described by a 
diffusion equation that governs the evolution of the probability distribution (p(X., x, t) of 
the random variable Xt (considered at point x in physical space): 

^ + u • V,^ = ^Axv^ - divx (V^uX (/p) + ^divx {F{X)ip) , (1.11) 

The expectation in (1.7) then becomes an average with respect to the probability measure 
(/7(X, x,t) dX: 

Tp{x, t) = ^(fx0 F(X) V9(X, X, t)dX - Id j . (1.12) 

We refer for example to [4, 8, 34] for more details on the physical background and more 
complicated models. 

A numerical simulation of the coupled system (1.5)^(1.8) is very expensive, since one 
needs to obtain the non-Newtonian stress tensor Tp at each space-time discretization node. 
Several approaches have been proposed in the literature [23, 28]. A first approach is a de- 
terministic micro-macro simulation. Here, one couples the Fokker-Planck equation (1.11)- 
(1.12) with the Navier-Stokes equations (1.5)-(1.6). The main drawback of these methods 
is their high computational cost, due to the high-dimensionality of the function (p (which 
depends on seven scalar variables (X.,x,t) in dimension 3). This difficulty becomes all the 
more severe when more refined models involving higher dimensional microscopic variables 
Xt are used to describe the polymers. Specialized techniques are currently being developed; 



see e.g. [1, 2, 7]. The micro-macro simulation can also be performed stochastically. One then 
discretizes the macroscopic fields (velocity, pressure, stress) on a mesh, and supplements 
the (macroscopic) discretization of the Navier-Stokes equations with a stochastic simulation 
of an ensemble of polymers using a discretization of the SDE (1.8), see [15, 27]. Methods 
have been proposed to obtain sufficiently low- variance results [6, 15, 20]. 

Due to the very high computational cost of micro-macro simulations, another route 
which has been followed (see e.g. [14, 16, 22, 32, 33, 35]) is to look for an approximate 
closure at the macroscopic level, namely a model of the form: 

dM 



dt 



u-V^M = H(M,V^u), (1.13) 

r(M), (1.14) 



which is close to the microscopic model (1.7)-"(1.8). Here M denotes an ensemble of macro- 
scopic state variables that depend on time and space. A basic example of such a macro- 
scopic model is the Oldroyd-B model [4], which is actually equivalent to the microscopic 
model (1.7)-(1.8) for a linear force F(X) = Fhook = X. In this case, one can obtain a 
closed equation on the so-called conformation tensor cr(t) = {aij{t))f-'^^, with d the num- 
ber of space dimensions, and (Jij{t) = {{Xi)t{Xj)t), in which {Xi)t, resp. {Xj)t, represent 
the corresponding component of X^. This yields the equation : 

dtTp + U • VxTp = VxUTp + TpVxU^ + — (Va;U + Va;U^) - ^Tp- 

On the other hand, for the FENE model, no equivalent closed macroscopic model is known, 
and one has to resort to approximate closures to obtain macroscopic equations (see Sec- 
tion 2.2). The basic idea is to approximate the polymer distribution by a so-called canonical 
distribution function, which is determined using only the macroscopic state variables M 
(typically low-order moments of the distribution). The microscopic evolution law (1.8) (or 
(1.11)) is then replaced by a set of equations (1.13) for the evolution of the macroscopic 
state variables M, combined with a constitutive equation (1.14) for the stress. While such 
approximate macroscopic models are desirable, at least from a computational point of view, 
it is however not always clear how to quantify the effects of the introduced approximations 
on the accuracy of the simulation, and how to choose the macroscopic state variables M. 

Recently, there has been quite some interest in the development of computational meth- 
ods that aim at accelerating micro-macro simulation using on-the-fly numerical closure 
approximations. We mention equation-free [24, 25] and heterogeneous multiscale methods 
(HMM) [10, 11]. In both approaches, a crucial step is to define an operator that generates a 
microscopic state corresponding to a given macroscopic state; this is actually equivalent to 
prescribing the closure approximation. This step is called lifting in the equation-free frame- 
work, and reconstruction in HMM. Inspired by these methods, the present paper studies 
in detail the question of lifting/reconstruction for the particular problem of micro-macro 



models for polymeric fluids; the procedm'e we propose, however, could be applied to many 
multiscale models. Specifically, we propose a computational procedure to reconstruct an 
ensemble of N polymers consistently with a given macroscopic state M, and we exam- 
ine the errors that are introduced in the macroscopic evolution by numerically enforcing 
closure upon the selected macroscopic state variables. For convenience of exposition and 
illustration, we restrict ourselves to one-dimensional simulations with pre-imposed (time- 
dependent) velocity fields, i.e. equations (1.7)-(1.8) with given u(x,i), at one specific point 
X in space. However, we emphasize that the numerical method can be used likewise for 2D or 
3D situations, as well as for the closure approximation for the coupled problem (1.5)-(1.8). 

The main contributions of the present paper are twofold: 

• From a modelling viewpoint, we propose a numerical closure strategy that enables to 
easily explore which sets of macroscopic state variables should be chosen to get good 
closure approximations. Various strategies are proposed and tested. 

• From a theoretical viewpoint, we show the relation between this numerical closure 
strategy and the so-called quasi-equilibrium method proposed in [19], which relies on 
an entropy minimization principle. 

The paper is organized as follows. In Section 2, we give some more detail on the FENE 
model and the existing literature on closure approximations. In Section 3, we propose 
a numerical closure approximation based on constrained SDE simulations [29], which is 
very fiexible, and enables to explore the error introduced by the closure for various sets 
of macroscopic state variables M. This numerical closure approximation is shown to be 
optimal in the sense that, when applied to a microscopic model which has an equivalent 
macroscopic model, it indeed yields the macroscopic model (Section 4). Moreover, we show 
in Section 5 that, in some specific cases, it is closely related to the closure approximation 
based on a quasi-equilibrium condition introduced in [19]. Finally, we test the numerical 
closure using a number of different strategies to define the macroscopic state variables 
M (Section 6). We first perform numerical experiments to assess the capability of the 
selected macroscopic state variables to recover the desired polymer distributions in strong 
flow regimes. Second, we study if the procedure is able to correctly capture macroscopic 
evolution. While accelerating microscopic simulation is not the primary purpose of the 
present paper, we give some remarks in this respect in Section 7, where we briefly discuss 
the main results and give some directions for future research. 



2 The FENE model and closure approximations 

2.1 FENE dumbbells: discretization and a one-dimensional version 

As mentioned above, we consider polymer simulations with FENE dumbbells subject to a 
pre-imposed (time-dependent) velocity field. Thus, in the remainder of the paper, unless 
explicitly stated otherwise, the force is the FENE force, see (1-10) : 



' FENE- 



Using the characteristic method to integrate the convective derivative in (1.8) (Lagrangian 
frame), the equations of interest reduce to: 



We 



(Xi®F(Xi))-Id , 



dXi 



n{t)y.t 



2We 



-F(X, 



f2.1) 



dt + 



/We 



dWi, 



where Xj now depends on the foot of the characteristic rather than on the Eulerian space 
position X, and k is the velocity gradient (along the trajectory). Unless stated otherwise, 
we will work with a one-dimensional version of this equation, 



4((X,F(A-,)>-l 



dX, 



K{t)Xt 



1 



2We 



F{Xt 



dt + 



/We 



dWu 



(2.2) 



keeping in mind that the algorithm, as well as its analysis and implementation extend 
straightforwardly to higher dimensions. Note that K{t) is here a given one-dimensional time- 
dependent function and F denotes a one-dimensional version of the FENE force, see (1.10), 
namely 

Such a one-dimensional framework has also been used in [22] for example, to assess the 
influence of the Peterlin approximation (see Section 2.2) on transient behaviour. 

Concerning discretization methods, we use a classical Euler-Maruyama scheme [26] with 
a Monte Carlo method: 



e 

We 



/ , N 

\ n=l 



1 



X 



n,k+l 



Xn,k^ 



K{t^)X 



n,k 



2We 



F{X 



n,k\ 



(2.3) 



6t + 



/We 



6t^ 



n,k 



where the indices n and k denote respectively realization index and time index, t^ = k6t 
and ^"' are i.i.d. normal random variables. 



For convenience, we introduce a short-hand notation for the discretization scheme of the 
SDE in (2.3), 

X'^^'=s^{X\>,{t'^),St), (2.4) 

where X = {X"^}^^-^ is the ensemble of A^ reaUzations, and K(i'^) indicates expUcitly the 
value of the velocity gradient in (2.3) that is considered over the time interval of size 5t. 

Theoretically, it can be shown that (for sufficiently large b), the norm of the end-to- 
end vector in (1.8) or (2.2) (recall that F = Ffene) cannot exceed the maximal value 
\/b [21]. However, the discretization scheme (2.3) might yield spring lengths beyond this 
maximal value. There are two ways to deal with this problem [34, Section 4.3.2]. The 
first is via an accept-reject method, in which, for each polymer, the state after each time 
step is rejected if |X|^ > (1 — v5t)&, and a new random number is tried until acceptance. 
Alternatively, one could use a semi-implicit predictor-corrector method. In this text, we 
choose the accept-reject strategy. 

2.2 Closure approximations for FENE dumbbells 

We now briefly discuss the derivation of closure approximations of the type (1.13)- (1.14) 
for the FENE model. 

One closure approximation is the Peterlin pre-averaging [5]. Here, one constructs an 
approximation for the FENE model by defining the spring force as (compare with (1.10)) 

V 

Ffene~p{X) = ^ _ {x'i)/b ' ^'^'^^ 

As a consequence, only the mean square length of the ensemble of polymers is constrained to 
remain smaller than v5, whereas the length of individual polymers may exceed this value. 
The interest of FENE-P dumbbells is that, as for Hookean dumbbells, a closed equation can 
be derived on the conformation tensor a = {X^), and thus a macroscopic model is obtained: 

dfcr + nVrO" = 2a'VrU —-—, — I . 

Wel-trCT/6 We , , 

/ X 2.6 



^ We Vl-tr(cj)/6 

It has been shown in [14, 22] that the Peterlin approximation has a profound impact on 
transient behaviour in complex flows, compared to the original FENE model. 

Let us now discuss more generally closure approximations of the type (1.13). For the 
sake of clarity, and without loss of generality, we restrict ourselves to the one-dimensional 
case (2.2). 

Consider starting from a number L of macroscopic state variables, M = {Mi}^^^, which 
are defined as configuration space averages of functions mi of the configuration Xt , 

Mi{t) = {miiXt)) . (2.7) 



The goal is to obtain a closed system of L evolution equations (1-13) for the state vari- 
ables M, complemented with a constitutive equation (1.14) for Tp as a function of these 
macroscopic state variables. 

Using Ito calculus, one can easily obtain the following equation of state for the macro- 
scopic state variables, 

™' .(<)^,?S^U.)V^^U.)^W)V^(|S™V (2.8) 



dt ^ '\ dX^ ' / 2We \ ' ' dX^ 7 2We\dX2 



MP MC MP 

in which the macroscopic state variables M^ ' ' account for hydrodynamic drag, connec- 
tor force and Brownian motion, respectively. Of course, in general, many of these macro- 
scopic state variables M^ ' ' are not functions of the initially chosen macroscopic state 
variables {Af;}^^^. One can write evolution equations for these new state variables, which 
in turn will create additional state variables but this procedure typically goes on endlessly. 
At some point, one has to stop, and try to approximate the state variables for which no 
evolution equation is available by writing them as a function of other (already available) 
state variables. By adding such closure relations, one obtains an explicit, but approximate, 
closed system of evolution equations. 

Any closed macroscopic model needs to (i) define the set of macroscopic state variables 
M = {Mi}i^^, and (ii) provide a way of evaluating the remaining state variables M^ ' ' 
in the evolution equation as a function of M. In the literature, item (i) is generally addressed 
by considering a hierarchy of even moments, i.e. Mi = {X ) where / = 1, . . . ,L (all the 
odd moments are zero for reasons of symmetry). Note that these become tensors in higher 
space dimensions. The corresponding evolution equations (2.8) are then given as: 

dMi , , , 1 r- ^(2/ - 1) , , 

with Mq = 1. In order to complete (ii), one needs to provide approximations for the 
new additional macroscopic state variables {Mp}._ . Note that, in particular, one of 
this new additional macroscopic variable M^ is also required to obtain the constitutive 
relation (1.14) for Tp. One strategy to approximate |Mp|._ is to propose a probability 
distribution (p-M{X) (called a canonical distribution function) that is parameterized by the 
selected macroscopic state variables, and to compute {Mp}, in the evolution equations 
(2.8) as the expectation with respect to this canonical distribution function. Note that 
'f'M.iX) depends on time only through the dependency of M on the time variable. The 
rationale behind this approach is that the better one can approximate the microscopic 
distribution function, the more reliable the obtained macroscopic model should be. 

In [32, 33], approximate closures for Mf are obtained by restricting the space of ad- 
missible distribution functions to linear combinations of L canonical basis functions. Based 
on this approach, several closures have been proposed; see [32] for more details on the 



one-dimensional setting (2.2) and [33] for the general three-dimensional case. A related 
approach is described in [9, 16, 36]. Another route is described in the following section. 

2.3 Quasi-equilibrium approximations 

A particularly interesting approach is proposed in [19]. It consists in defining a so-called 
quasi-equilibrium canonical distribution function Lp^ via a constrained entropy optimiza- 
tion problem: 

^2a = argmin hpln ( ) , (2.10) 

where Vt-^ is defined as the set of all probability density functions, for which the average of 
mi is indeed M/: 



n 



M 



|v9(X),(^>0, fipiX)dX = l, fmiiXMX)dX = Mi,l = l,...,L\. (2.11) 



In (2.10), ifeq is defined as the equilibrium distribution for the polymer configuration, for 
zero velocity field. In particular, for FENE dumbbells, it writes: 

^,,{x) = z-'{i-x'/bf\ 

where Z= I {l - X^ /b)^^^ dX. 

J\X\<Vb 

The rationale behind this approximation is to assume a separation of time scales be- 
tween the (supposedly fast) relaxation towards the quasi-equilibrium distribution and the 
(supposedly much slower) evolution of the macroscopic state variables. 

An explicit expression of the solution to (2.10) can be obtained as: 

<f'^{X) = Z^^,,{X)e^p\S2Ximi{X)j , (2.12) 

where Z^ = / ipeq{X) ex.p N^A/r?T,i(X) dX and the set of Lagrange multipliers A = 

{Xi}ILi are determined by the constraints / mi{X)ip^ {X) dX = Mi. 

While the Lagrange multipliers depend only on the macroscopic state M, the relation 
A(M) can often not be obtained analytically. Therefore, in [19], a numerical procedure is 
proposed to simulate the resulting closed macroscopic model. We will show below (see Sec- 
tion 5) that the numerical closure approximation technique that we propose (see Section 3) 
is closely related to this method, and that it may be considered (for a slightly modified ver- 
sion) as a different numerical strategy to obtain quasi-equilibrium closure approximations. 



3 Numerical method 

In this section, we propose to mimic the evolution of the corresponding unavailable macro- 
scopic model via a coarse time-stepper [24, 25]. 

3.1 The lifting and restriction operators 

Consider the evolution of an ensemble of polymers in a pre-imposed velocity field and 
define a set of macroscopic state variables M which are believed to represent the underlying 
(microscopic) polymer distribution sufficiently accurately. We introduce two operators that 
make the transition between microscopic and macroscopic state variables. We define a 
lifting operator, 

c-.m^a:, (3.1) 

which maps a macroscopic state to an ensemble of N polymer configurations, and the 
associated restriction operator, 

n-.X^M, (3.2) 

which maps an ensemble of configurations to the corresponding macroscopic state. Note 
that we directly define the method at the discrete level over an ensemble of A^ configurations 
(after Monte Carlo discretization). For a discussion in the limit of an infinitely large number 
of polymer configurations, we refer to Section 3.3. 

The restriction operator is readily defined using an empirical mean: 

1 ^ 
n{X) = {Ml = ni{M)}{L^ with TZiiX) = ^ J^rn^X-) for / = 1, . . . ,L, (3.3) 

n=l 

where, we recall, X = {X^}^^-^ denotes the ensemble of configurations. 

In the lifting step, we need to sample a reconstructed polymer distribution function, 
consistently with the given macroscopic state M(t*) obtained at time t* . To this end, 
we perform a constrained simulation of an ensemble of polymers until equilibrium, subject 
to the constraint that the macroscopic state remains constant and equal to M(t*). More 
precisely, the constrained algorithm writes [29]: 

L 

ti (3.4) 

with A G M^ such that 7^^(A'"'+^) = Mi{t*) for / = 1, . . . , L. 

It thus consists successively in an unconstrained Euler-Maruyama step, followed by a pro- 
jection step to satisfy the constraint. In each constrained time step, the projection is done 
by solving the nonlinear system 

lZi{X'^+^{k;X'^,6t))=Mi{t*), for / = !,..., L, (3.5) 
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for the unknown Lagrange multipliers A using Newton's method. In (3.5), we have made 
explicit that the state Af™"*"^ depends on the unknown Lagrange multipliers, as well as on 
(known) X"^ and 5t. During the constrained simulation, an accept-reject strategy is applied 
on the combined evolution and projection operation, i.e. if, during projection, the state of 
a polymer would become unphysical, we reject the trial move in the unconstrained Euler- 
Maruyama step and repeat the time step for this polymer, after which the projection of the 
ensemble is tried again. 

The lifting operator is then defined as the ensemble X^°° for a sufficiently large time 
index moo, which is chosen such that (3.4) has reached an equilibrium distribution, 

£(M) = A'™-. (3.6) 

We will detail further on how moo is determined numerically when describing the computa- 
tional experiments. For a precise definition of the lifting operator in terms of distributions 
(in the limit of an infinite number of configurations), we refer to Section 3.3. 

Of course, by construction one has the consistency property 

7^ o £ = Id. 

3.2 The numerical closure algorithm 

Let us now make precise the complete algorithm. Given an initial condition for the macro- 
scopic state variables M(t*) at time t*, one time step of the coarse time-stepper consists of 
a three-step procedure: 

(i) Lifting, i. e. the creation of initial conditions 

X{t*) = £(M(t*)) 
for the microscopic model, consistently with the macroscopic state M(i*) at t* . 

(ii) Simulation using the microscopic model over a time interval \t* , t* + K5t\ , where K 
is the number of time steps, to get X{t* + K5t): for A: = 0, . . . , ii' — 1, 

X{t* + {k + l)5t) = sx {X{t* + k6t), K{t* + k6t),6t) . 

(iii) Restriction, i.e. the observation (estimation) of the macroscopic state at t* + K5t: 

M{t* + K6t) = n{X{t* + K5t)). 

In the following, we denote 

At = K5t. 



11 



During the restriction step, the ensemble X(t* + At) is also used to get an estimate of 
the new value of the stress 

r,{f + At) = ^ f 1 f; X^{t* + At) F{X-{t* + At)) - 1 j . 

3.3 The lifting and restriction operator in the continuous Hmit 

The lifting and restriction operators which have been defined above depend on three dis- 
cretization parameters: A^ which is related to the Monte Carlo discretization (the operators 
have been defined for a finite ensemble of configurations), 6t which is related to the time 
discretization in (3.4), and rrioo which should be sufficiently large to reach a stationary state 
in (3.4). In this section, we introduce the limiting operators C and TZ obtained in the limit 
N — )• oo, (^t — )• and niooSt — )• oo. 

Note first that these operators are well-defined in terms of the probability distribution 
ip, rather than ensembles of configurations. More precisely, the lifting operator C con- 
sists in constructing a probability distribution ^^ consistently with the macroscopic state 
variables M (using the notation of Sections 2.2-2.3), 

£(M) = ^^^(X), (3.7) 

in which the superscript NC stands for numerical closure. Likewise, the restriction operator 
7^ reduces a distribution to macroscopic state variables. 

The restriction operator TZ is simply an averaging operator, which computes the averages 
of rui with respect to the distribution (p (compare with (3.3)): 

TZiip) = {Ml = niiip)}iLi with TZii^) = fmiipioi:l = l,...,L, (3.8) 

On the other hand, the lifting operator C is more involved to define. When considering 
the continuous-in-time version of (3.4) in the limit of an infinite number of configurations, 
A^ — 7- oo, it can be seen to be given by the one-dimensional marginal of the stationary state 
of the associated Fokker-Planck equation. 

Let us make this statement precise. For a fixed value A^, the numerical scheme (2.3) 
is a discretization of the following constrained Stratonovitch SDE on the ensemble Xt = 
{X^j^^i (see [29] and [30, Chapter 3]): 



dXt = P(A'i) 



K(t*)Xt —Fi^t) 

^ ' 2We ^ ' 



dt + — Lp(A't) o dWt, (3.9) 

vWe 



where, with a slight abuse of notation, F{Xt) = (P(^t"))„=i, and Wt represents an A^- 
dimensional Brownian motion. The projection operator P{Xt) is defined by: 

L 

P{X) = Id - J^ Grj{X)Vxni{X) €5 Vxn,{X) 
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with G- ■ iX) the inverse of the Gram matrix: 

G^,J{X) = VxMX) ■ VxTljiX) 

and o denotes the Stratonovitch product. If we denote 

S(M) = {X, n{X) = M} (3.10) 

the submanifold of X at fixed values of the macroscopic state variables, then P{X) is the 
orthogonal projection operator onto the tangent space T:f S(M) of S(M) at point X. Thus, 
if Xo G S(M), then, for all t > 0, Aft G S(M). 

Let us denote ip {t,dX) the distribution of Xt satisfying (3.9). Note that the compo- 
nents of Xt have all the same law, for symmetry reasons. Let us introduce the marginal of 
■0^ in the first variable: 

ij^{t,X^)dX^= [ i;^{t,dX\...,dX^). (3.11) 

Jx^,...,xN 

Then, ip^ is defined as: 

</'m''(^) = lim lim <(t, X). (3.12) 

By a law of large numbers, it is expected that this distribution ^^ is consistent with the 
fixed values of macroscopic state variables M: 

where Qm is defined by (2.11). 

We will discuss in Section 5 how to get an analytical expression for v'm ' ^^ least in 
some specific cases. 

3.4 Choice of the macroscopic state variables 

For the FENE model, it appears that the first even moment (X^) is not sufficient to char- 
acterize the polymer distribution, and additional macroscopic state variables are needed. 
We will consider the macroscopic level to be determined by L macroscopic state variables, 
M = {Mi}JL^, and we consider the following strategies to select Mi, I = 1, . . . ,L. 

Strategy 1. We consider a hierarchy of even moments of increasing order, 

Mi = {X^'), l = l,...,L. (3.13) 
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Strategy 2. We consider a hierarchy of even moments of increasing order, and supplement 
the set of macroscopic state variables with the additional moments that appear in the 
corresponding evolution equations (2.9), 

\ M^/2+« = Mf = 2l{F{Xt) Xf-'), 

for 1 < / < L/2 where L is assumed to be even. For FENE dumbbells, it can easily be 
checked that 

and that all Mp, I > 1 can be written as linear combinations of Mi, 1 = 1,..., L/2 and Tp. 
Hence, this choice is equivalent to taking 

'Mi = {Xf), l = l,...L-l 

f f // X2 \ \ (3.15) 

A4=r, = -((A-,FW))-l) = -((^-^\ '^ 



where L = L/2 + 1 denotes the number of linearly independent macroscopic state variables. 

Strategy 3. We again start from Mi = (Xf). To add state variables, we write down the 
evolution equation for Mi, i.e. (2.9) with / = 1, and add all macroscopic state variables 
that appear in this equation. In this case, this amounts to adding the variable M2 = Mf . 
We continue by writing down the evolution equation (2.8) for M2, which, in turn, reveals 
additional state variables M2 ' ' . Some elementary algebra shows that we obtain four 
linearly independent macroscopic state variables: 

Mi = (x2), M2 = ^ ^_^,^J -1, (3.16) 

as above, and additionally 

Note that these same macroscopic state variables would also show up after simplification 
by applying this procedure starting from the choice Mi = Tp. If additional moments are 
desired, one could continue by writing down evolution equations for M3 and M4 and add 
the moments that appear in those equations, but we will not consider that in the remainder 
of the text. 



4 A consistency result for FENE-P dumbbells 

To check the consistency of the whole procedure, let us apply the numerical closure ap- 
proximation to the case of FENE-P dumbbells (namely using the spring force (2.5)). In 
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this case, it is known that there exists a macroscopic equivalent model and the question is 
thus: do we recover this macroscopic model using the numerical closure procedure ? We 
first derive a theoretical result, which we subsequently illustrate numerically. 

4.1 A simple remark 

Let us consider the FENE-P model, with the above numerical closure approximation method 

applied using only one macroscopic state variable M = {Xf). Note that the stress Tp is 

defined in terms of M as 

€ f M 



^ We V 1 - M/b 

As mentioned above (see (2.6)), for the microscopic model (2.2), M satisfies a closed 
equation: 

dtM = 2KM-— ^^ + —. (4.1) 

Wel-M/6 We ^ ^ 

We now make a simple observation to show that the numerical closure approximation 
(in the limit of zero discretization errors) reproduces this macroscopic dynamics. We refer 
to the notation of Section 3.2. For a given value of M{t*) at time t* , the lifting step (i) 
creates an ensemble of configurations with, by construction, a law y^^m*) — C{M{t*)) such 

that / X'^ip^,f^^,JX) dX = Af (t*). But then, the simulation step (ii) will indeed propagate 
M according to (4.1) (which is deduced from (2.2) by a simple Ito calculus). Thus, after 
the restriction step (iii), the correct values for M are recovered. 

In conclusion, if there exists a closed macroscopic equation for the stress, the proposed 
numerical closure approximation indeed recovers this macroscopic evolution as soon as the 
appropriate macroscopic state variables are selected. 

4.2 Numerical illustration 

We consider one-dimensional FENE-P dumbbells, governed by (2.2), in which the spring 
force F{X) = Fpene-p{X) is given by (2.5) with nondimensional parameters b = 49, 
We = 1 and e = 1. As in [22], we prescribe the velocity field 

K{t) = 100 t (1 - t) exp(-4t). (4.2) 

The microscopic model (2.2) is discretized via the Euler-Maruyama method with time step 
5t = 10-2. 

4.2.1 Lifting 

To illustrate that the macroscopic variable M = {Xf) uniquely determines the polymer 
distribution, we perform the following experiment. We first simulate an ensemble of A^ = 10^ 
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Figure 1: Polymer distribution for FENE-P dumbbells during constrained simulation. Shown are 
the polymer distribution before the restriction at i = 0.3 (the reference distribution), and at several 
time instances during a constrained simulation starting from a uniform initial distribution. (The 
non-uniform appearance of the initial condition is due to artifacts of the binning.) Parameters of 
the simulation are given in the text. 

FENE-P dumbbells, subject to the velocity gradient K{t) over the time interval t G [0,0.3]. 
As the initial condition, we take the equilibrium polymer distribution in the absence of 
flow. At t = 0.3, we obtain M* = M{t = 0.3) via restriction; the corresponding polymer 
distribution is kept as the reference distribution. Next, we initialize a new ensemble of 
polymers consistently with the macroscopic state M* using a uniform distribution. We then 
perform a constrained simulation (3.4) using the same time-step 5t over the constrained time 
interval [O,moo'^i] = [0, 50]. The results are shown in Figure 1. 

We see that the distribution of the constrained simulation converges towards the dis- 
tribution of the original simulation, indicating that the first even moment M is indeed 
sufficient to represent the original polymer distribution, and also that the constrained sim- 
ulation recovers this distribution. 

Note, however, that this experiment reveals an important property of FENE-P dumb- 
bells. While the manifold consisting of Gaussian distributions with zero mean is invariant, 
there is no strong time-scale separation between the relaxation of arbitrary distributions 
with given second moment towards the Gaussian distribution and evolution of this second 
moment itself. This can be concluded by noting that one needs to simulate the constrained 
SDE over a time interval of length 50 to reach the stationary distribution, whereas the 
macroscopic state variable evolves significantly on considerably shorter time-scales, see also 
the next experiment. This was also observed in [17]. 
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Figure 2: Evolution of the first even moment M and stress Tp for an ensemble of FENE-P dumbbells 
during complex flow. Left: (Af, Tp) phase plane view. Right: temporal evolution. Shown are a full 
microscopic simulation (reference), and simulations using a coarse time-stepper for different values 
of the macroscopic time-step. Simulation parameters are given in the text. 

4.2.2 Coarse time-stepping 

We now look into the evolution of the numerical closure with respect to the full microscopic 
simulation. To this end, we simulate an ensemble of A^ = 2-10^ FENE-P dumbbells, starting 
from the equilibrium distribution feq in the absence of flow, up to time t = 2. All numerical 
parameters are the same as above. In particular, K{t) is again given by (4.2). We compare 
this reference simulation with a number of simulations using the coarse time-stepper with 
different values of the time step At = K5t. In this experiment, the lifting step amounts to 
freezing physical time and performing a constrained simulation that is consistent with M. 
The constrained simulations are performed over a time interval of size lOOAt. The results 
are shown in Figure 2. We see that the results are nearly identical for all values of At and 
the results nearly coincide with the reference simulation. This is to be expected. Indeed, 
since M completely determines the polymer distribution, a simulation constrained upon M 
will not alter this distribution, see Section 4.1. 

5 Comparison of numerical closure with quasi-equilibrium 
method 

In this section, we compare the proposed numerical closure approximation (described in 
Section 3) with the quasi-equilibrium method proposed in [19] (described in Section 2.3). 
In particular, we show that the quasi-equilibrium method, as proposed in [19], is equiv- 
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alent to the numerical closure approximation, when the velocity gradient n{t*) is taken 
to zero in (3.4). To prove this result, we need to show that the canonical distribution 
'^M reconstructed from the quasi-equilibrium method (see Equation (2.12)) is the same as 
the distribution v'm^ reconstructed from the lifting procedure through the operator C (see 
Equations (3.7) and (3.12)). 

Let us consider the microscopic model (1.7)-(1.8), with a general force F which derives 
from a potential IT: 

F = vn, 

so that the equilibrium distribution (for zero velocity field) is 

ipeq = Z"^exp(-n), 

where Z = Jexp(— H). Let us consider a fixed given set of macroscopic state variables 
M, and, for the sake of simplicity, let us assume that L = 1 (only one macroscopic state 
variable M is considered). 

From the quasi-equilibrium method, the reconstructed distribution is (see Equation (2.12)): 
ip'^fiX) = Z2f exp (-n(X) + Xm{X)) , (5.1) 

where Z^ = / exp (— n(X) + Xm{X)) dX and the single Lagrange multiplier A is deter- 
mined by the constraint / m{X)ip^,j (X) dX = M. 

Let us now consider the numerical closure approximation described in Section 3, with 
K,{t*) = in (3.4). In this case, since K{t*) = in (3.9), the stationary distribution for (3.9) 
has a simple expression: 

N 

^^(oo,dA') = (Z^)-i J] exp(-n(X"))das(M), 

n=l 

where o"2(m) is the Lebesgue measure on the submanifold S(M) defined by (3.10). We 
refer for example to [29] or [30, Proposition 3.20]. Then, the marginal ip^ {oo, X) is defined 
through (see (3.11)): 

ij^{oo,X^)dX^= [ tP^{oo,dX\...,dX^), (5.2) 

Jx'2,...,XN 
and the reconstructed distribution from the numerical closure approximation is (see (3.12)): 

(^j^^(X)= hm <(oo,X). (5.3) 

The main mathematical result of this work is the following: 

Proposition 5.1. The reconstructed distributions obtained through the quasi-equilibrium 
method, and the numerical closure approximation method with zero gradient velocity field 

are the same: 

^QE NC 

^M = Vm ■ 
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This proposition is a corollary of a general result about the equivalence (for an infinite 
number of particles) of the canonical ensemble and the microcanonical ensemble in statistical 
physics. We cite a result from [3, Theorem A. 5. 5], see also [13, Theorem 3.4]: 

Theorem 5.2. Let a be a probability measure on M'^ and let us consider Y^, . . . ,y^ i.i.d. 
random variables with law a, and introduce a function g : M — )■ M. Let us now define two 
probability measures: 

• The conditional measure 

( ' ^ 



V Y. i^y" 



N 

n=l 



of the vector {Y^, ..., Y ) conditionally to — X]n=i qO^^) — ^■ 



• The probability measure 

ax{dy) = Z~^ exp(Ag(y)) a{dy), 

where Z\ = J exp{Xq{y)) a{dy). Let us assume that A and z are related through the relation: 

q{y)a\{dy) = z. 
Then, one has: for any test function F : M'^ — )• M, 

hm f F{y^)vndy\...,dy'')= I F{y^)a^{dy^). 



N- 



To apply Theorem 5.2 to prove Proposition 5.1, we set a to be the equilibrium distri- 
bution ifeq, q = rn, and z = M. Then ax = ipj^j , and it remains to show that 

,l^,„JVC'^„,l^ T ,1 



lun^l Fiy')u(: {dy\...,dy'^) = J Fiy')^Z''iy') 
This is stated in the following lemma: 



Lemma 5.3. Let us consider the notation of Theorem 5.2 and assume that the measure a 
has a density a: 

a{dy) = a{y) dy. 

Let us introduce the probability measure 

.N\ 



'{z) 



where S^(-z) = {{y^, ■ ■ ■ ,y^), jf Sn=i liv"') — ^} '^'^d (Jj^nu) ^^ ^^e Lebesgue measure on 



the submanifold T,^ (z) . Then, 

u^^^^^idy\...,dy'') = \\VQ''\\u('^{dy\...,dy''), (5.4) 

where Q^ {y^ , • • • , y^) = -^ Ylin=i Qiv"')- Moreover, 

hm f F{y^)u('^{dy\...,dy'')= Yira f F{y^)u^^,Jdy\ . . . ^dy"") . (5.5) 
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Proof. The proof of (5.4) is based on the co-area formula, see for example [30, Eq. (3.14)]. 
Then, to prove (5.5), one notice that, if Y^, . . .Y denotes random variables distributed 
according to the conditional probability measure z/,^, one has: 

j,( U N f,i , 7V^ {F{Y^)\\VQ'^\\{Y\...,Y^)) 



^Etillv<?P(y«; 



By a law of large numbers (see for example [3, Theorem A. 5.4] or [13, Theorem 3.5]), 

1 ^ f 

— N^ ||V(7|| (y") converges in probability to / ||Vg|| doA) and thus, Slutsky lemma enables 

n=l -^ 

to conclude. D 

This concludes the proof of Proposition 5.1, since with the notation introduced above 
{a{dy) = ipeqiy)dy, q = m, and z = M) 

i^^N^,)idy\ ..., dy^) = ^'^(oo, dy\..., dy^). 

A few remarks are in order. First, in dimension 1, the fact that the drift in the SDE 
derives from a potential is not a restrictive assumption, so that the quasi-equilibrium pro- 
cedure could also be applied when taking into account a non-zero K{t*). However, this 
assumption is indeed restrictive in dimension greater than one: for non-symmetric ^(t*), 
the drift in (2.1) is not the gradient of a potential. In this case, the numerical closure 
approximation procedure still applies, but it is unclear how it would be related to a quasi- 
equilibrium method. In some sense, the numerical closure method can thus be seen as a 
generalization of the quasi-equilibrium method, which takes into account the velocity gradi- 
ent in the lifting procedure. In fact, the numerical closure procedure can be seen as a simple 
alternative to simulate the quasi-equilibrium closures that, unlike the numerical procedure 
in [19], does not require transformations from moments to Lagrange multipliers and vice 
versa, which might be difficult to perform. 



6 Numerical illustrations for FENE dumbbells 

In this section, we perform some numerical experiments to explore the behaviour of the 
numerical closure procedure using the strategies for macroscopic state variable detection 
that were outlined in Section 3.4. 
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Figure 3: Lifted polymer distributions for FENE dumbbells as a function of the number of macro- 
scopic state variables using strategy 1. We plot a reference polymer distribution, that is obtained 
by microscopic simulation up to time t* , as well as the equilibrium polymer distributions after con- 
strained simulation using L — 1, ... ,4 even moments. Shown are the results for t* = 0.5 (top left), 
t* — 1 (top right), t* = 1.5 (bottom left) and t* — 2 (bottom right). Simulation parameters are 
given in the text. 

6.1 Strategy 1: Even moments as macroscopic state variables 

6.1.1 Lifted configuration distributions 

We simulate an ensemble of A^ = 5 • 10^ FENE dumbbells, subject to a constant velocity 
gradient K{t) = 2 over the time interval t £ [0,t*], with t* = 0.5,1,1.5,2 (startup of 
"elongational" flow). We use nondimensional parameters 5 = 49 and We = e = 1, and 
choose 6t = 2 ■ 10~^. As the initial condition, we take the equilibrium polymer distribution 
in the absence of flow. As the macroscopic state variables, we take the first L even moments. 
At t = t*, we obtain M* = TZ{X*) via restriction; the corresponding polymer distribution 
is kept as the reference distribution. Starting from X* , we then perform a constrained 
simulation under the constraint that TZ{X) = M*, using the same time-step 6t, until the 
polymer distribution equilibrates. Figure 3 shows the constrained equilibrium polymer 
distributions for a range of values of L. We see that, as the number of macroscopic state 
variables increases, the difference decreases between the constrained equilibrium distribution 
and the reference distribution, indicating that this distribution is captured more accurately 
when more macroscopic state variables are used. 
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Figure 4: Evolution of the stress tensor Tp throughout constrained simulation using strategy 1 with 
L ~ 1 (top left), L — 2 (top right), L — 3 (bottom left) and L = 4 (bottom right). Solid lines are 
obtained using the procedure outlined in Section 3; dashed lines correspond to the quasi-equilibrium 
approximation. Simulation parameters are given in the text. 

6.1.2 Relaxation to equilibrium and comparison with quasi-equilibrium ap- 
proach 

We now repeat the above experiment with A^ = 2000 particles and t* = 1, and plot the 
evolution of the polymer stress Tp as a function of time. All other simulation parameters are 
as above. Moreover, to obtain the corresponding result for the quasi-equilibrium method of 
[19], we perform the same experiment, but now with «;(t) = throughout the constrained 
simulations. We ensured that both constrained simulations were performed using the same 
random numbers. The results are shown in Figure 4. The figures clearly show a relaxation 
towards the stress value that corresponds to the lifted polymer distribution. This fact 
can be used to detect when the constrained simulation has equilibrated, and hence to 
determine the parameter moo that was introduced when defining the lifting operator in 
Section 3. When using the other strategies to determine the hierarchy of macroscopic state 
variables, Tp belongs to the set of macroscopic state variables, and therefore does not change 
during relaxation. However, in similar experiments, not reported here, we observed similar 
behaviour when monitoring the first even moment that was not constrained. 

Moreover, when the number of macroscopic state variables increases, the stress Tp that 
corresponds to the lifted distribution approaches the stress associated with the distribution 
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that corresponds to the initial condition of the constrained simulation. This observation is in 
agreement with the previous experiment, where we showed that the distributions themselves 
approach the initial distribution of the constrained simulation when more moments are 
taken into account. Hence, monitoring the evolution of Tp during constrained simulation 
can be used to determine whether the currently used set of macroscopic state variables is 
sufficient. Finally, concerning the relation between the numerical closure and the quasi- 
equilibrium approximation, we see that the difference between the two approaches is not 
really large; however, this difference remains of the same order of magnitude, independently 
of the number of macroscopic state variables included. 

6.1.3 Coarse time-stepping 

We now look into the evolution of the numerical closure with respect to the full microscopic 
simulation, again using K{t) = 2. To this end, we simulate an ensemble of A^ = 2000 
FENE dumbbells, starting from the equilibrium distribution in the absence of flow, up to 
time t = 4. All parameters are the same as above. We compare this reference simulation 
with a simulation via the coarse time-stepper, using a range of values for the number L 
of macroscopic state variables; here, the macroscopic time-step is equal to one microscopic 
step 6t, i.e. K = 1. In this experiment, the lifting step amounts to freezing physical 
time and performing a constrained simulation that is consistent with M. The constrained 
simulations are performed until equilibrium of the distribution is reached (here using ttt-oo = 
50 constrained time steps of size 5t)] all simulations were verified to have converged with 
respect to the number of constrained time-steps. The results are shown in Figure 5. We 
clearly see that the approximation improves as a function of the number of moments that 
are included at the macroscopic level. Other experiments, not reported here, indicate that 
the higher K,{t), the higher the number of macroscopic state variables that needs to be 
considered. These results are in line with the conclusions in [19], where analytical (quasi- 
equilibrium) closures were obtained via an entropy maximization principle. 

Finally, we consider an ensemble of A^ = 2000 FENE dumbbells subject to the time- 
dependent flow field 4.2, and again look at a coarse time-stepper in which the macroscopic 
state is represented with an increasing number of even moments. For this test, moo = 100; 
all remaining simulation parameters are as above. The results are shown in Figure 6. The 
conclusions for this experiment are similar. Note that a macroscopic description with only 
one moment cannot capture the hysteretic effect of the FENE dumbbells. 

6.2 Strategy 2: Adding the stress tensor as a macroscopic variable 

One particular advantage of the numerical closure strategy described here is that one can 
readily consider the effect of considering more complicated moments in the set of macro- 
scopic state variables. In this section, we repeat the above experiments, now considering the 
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Figure 5: Evolution of first even moment Mi (left) and stress r^ (right) for an ensemble of FENE 
dumbbells during startup of elongational flow. Shown are a full microscopic simulation (reference), 
and simulations using a coarse time-stepper for different numbers L macroscopic state variables 
using strategy 1. Simulation parameters are given in the text. 






Figure 6: Evolution of first even moment Mi and stress Tp for an ensemble of FENE dumbbells 
during complex flow. Left: (Mi,Tp) phase plane view. Right: temporal evolution. Shown are a full 
microscopic simulation (reference), and simulations using a coarse time-stepper for different numbers 
of macroscopic state variables using strategy 1. Simulation parameters are given in the text. 
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Figure 7: Lifted polymer distributions as a function of the number of macroscopic state variables. We 
plot a reference distribution, i.e., the polymer distribution after a microscopic simulation up to time 
t* , as well as the equilibrium polymer distributions after constrained simulation with L = 2, . . . , 5 
moments using strategy 2. Shown are the results for t* = 0.5 (top left), i* = 1 (top right), t* = 1.5 
(bottom left) and t* — 2 (bottom right). Simulation parameters are given in the text. 

first L — 1 even moments, supplemented with the stress Tp itself as a macroscopic variable, 
i.e., M = iMi)JL^ with Mi = (X^') for 1 < / < L - 1, as before, and Ml = Tp. 

6.3 Lifted configuration distributions 

We again simulate an ensemble of A^ = 5 • 10^ FENE dumbbells, subject to a constant 
velocity gradient K,{t) = 2 over the time interval t G [0,t*], with t* = 0.5,1,1.5,2 (startup 
of elongational flow) and obtain M* = TZ{X*) via restriction; the corresponding polymer 
distribution is kept as a reference distribution. We perform a constrained simulation, start- 
ing from X* , under the constraint that TZ{X) = M* using the same time-step 6t, until the 
polymer distribution equilibrates. Figure 7 shows the constrained equilibrium polymer dis- 
tributions for a range of values of L. Compared to the case when only even moments were 
used, we see that adding Tp as a macroscopic variable dramatically improves the obtained 
equilibrium distributions, and less moments may suffice to characterize the distributions. 
However, when L = 2 and L = 3, we see a peculiar artifact in the distributions, in the sense 
that we obtain an increase of the number of polymers with near-maximal length (a small 
second peak in the distributions on the right). This results in high probability of rejections 
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Figure 8: Evolution of the first even moment Mi (left) and stress Tp (right) for an ensemble of FENE 
dumbbells during startup of elongational flow. Shown are a full microscopic simulation (reference), 
and simulations using a coarse time-stepper for different numbers of macroscopic state variables 
using strategy 2. Simulation parameters are given in the text. 

throughout the constrained simulation. 

6.3.1 Coarse time-stepping 

We now look at the evolution of the numerical closure with respect to the full microscopic 
simulation, again using K(t) = 2. We simulate an ensemble of A^ = 2000 FENE dumbbells, 
starting from the equilibrium distribution in the absence of flow, up to time t = 4 and com- 
pare this reference simulation with a number of simulations using the coarse time-stepper 
with a different number p macroscopic state variables (L — 1 even moments, supplemented 
with the stress tensor Tp). As before, we choose the macroscopic time-step equal to one 
microscopic step 5t, i.e., K = 1; all other parameters are also chosen as above. We lift by 
freezing physical time and performing a constrained simulation that is consistent with M 
until equilibrium of the distribution is reached (here using moo = 50 constrained time-steps 
of size 5t) ; all simulations were verified to have converged with respect to the number of con- 
strained time-steps. The results are shown in Figure 5. Also here, we see an improvement; 
the result of the complex flow experiment is shown in figure 9. 

6.4 Strategy 3: Cascading from the equation of state for Tp 

Finally, we repeat the above experiments, now considering the moments to be determined 
by Strategy 3. 
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Figure 9: Evolution of first even moment Mi and stress Tp for an ensemble of FENE dumbbells 
during complex flow. Left: (Afi,Tp) phase plane view. Right: temporal evolution. Shown are a full 
microscopic simulation (reference), and simulations using a coarse time-stepper for different numbers 
of macroscopic state variables using strategy 2. Simulation parameters are given in the text. 

6.4.1 Lifted configuration distributions 

We again simulate an ensemble of A^ = 5 • 10^ FENE dumbbells, subject to a constant 
velocity gradient K(t) = 2 over the time interval t € [0,t*], with t* = 0.5,1,1.5,2 (startup 
of elongational flow) and obtain M* = TZ{X*) via restriction; the corresponding polymer 
distribution is taken as the reference distribution. We perform a constrained simulation, 
starting from X* , under the constraint that TZ{X) = M* using the same time-step 6t, 
until the polymer distribution equilibrates. Figure 10 shows the constrained equilibrium 
polymer distributions for an increasing number of macroscopic state variables. Compared 
to the previous two strategies, we here observe very good agreement with the reference 
distribution with less macroscopic state variables. 

6.4.2 Coarse time-stepping 

We now look at the evolution of the numerical closure with respect to the full microscopic 
simulation, again using K{t) = 2. We simulate an ensemble of A^ = 2000 FENE dumbbells, 
starting from the equilibrium distribution in the absence of flow up to time t = 4 and 
compare this reference simulation with a number of simulations using the coarse time- 
stepper with a different number L macroscopic state variables as above. As before, we 
choose the macroscopic time-step equal to one microscopic step 5t, i.e., K = 1; all other 
parameters are also chosen as above. We lift by freezing physical time and performing 
a constrained simulation that is consistent with M until equilibrium of the distribution is 
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Figure 10: Lifted polymer distributions as a function of the number of macroscopic state variables 
using strategy 3. We plot a reference polymer distribution, i.e., the polymer distribution after a mi- 
croscopic simulation up to time t* , as well as the equilibrium polymer distributions after constrained 
simulation using the indicated macroscopic state variables. Shown arc the results for t* = 0.5 (top 
left), t* — 1 (top right), t* = 1.5 (bottom left) and t* = 2 (bottom right). Simulation parameters 
are given in the text. 
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Figure 11: Evolution of the first even moment Afi (left) and stress Tp (right) for an ensemble 
of FENE dumbbells during startup of elongational flow. Shown are a full microscopic simulation 
(reference), and simulations using a coarse time-stepper for different numbers of macroscopic state 
variables using strategy 3. Simulation parameters are given in the text. 

reached (here using moo = 50 constrained time-steps of size 5t) ; all simulations were verified 
to have converged with respect to the number of constrained time-steps. The results are 
shown in Figure 11. Also here, we see the improvement; the result of the complex flow 
experiment is shown in Figure 12. 

7 Conclusions and discussion 

We proposed a numerical closure strategy that enables to easily explore which sets of macro- 
scopic state variables should be chosen to get good closure approximations for the kinetic 
simulation of polymeric fluids. The method involves the reconstruction of a polymer distri- 
bution as the constrained equilibrium of a microscopic Monte Carlo simulation, constrained 
upon the desired macroscopic state. The resulting algorithm is very flexible, and enables to 
explore the error introduced by the closure for various sets of macroscopic state variables 
M. We showed that this numerical closure approximation is optimal, in the sense that, 
when applied to a microscopic model which has an equivalent macroscopic model, it indeed 
yields the macroscopic model. Moreover, in some specific cases, the approach is shown to be 
closely related to the closure approximation based on a quasi-equilibrium condition. While 
the exposition in the present paper was restricted to the one-dimensional case, extensions 
to higher space dimensions are straightforward. 

The procedure straightforwardly enables to test hypotheses on which macroscopic state 
variables should be included to build good closures. We have examined three strategies 
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Figure 12: Evolution of first even moment Mi and stress Tp for an ensemble of FENE dumbbells 
during complex flow. Left: (Afi,Tp) phase plane view. Right: temporal evolution. Shown are a full 
microscopic simulation (reference), and simulations using a coarse time-stepper for different numbers 
of macroscopic state variables using strategy 3. Simulation parameters are given in the text. 

to define a hierarchy of macroscopic state variables. Our numerical experiments indicate 
that, at least for the cases considered in this paper, fewer macroscopic state variables are 
required to obtain accurate results when choosing a strategy that adds macroscopic state 
variables based on the unknowns that appear on the right-hand side of an Ito calculation 
for the already included state variables (Strategy 3 in this text). Moreover, the experiments 
in section 6.1.2 indicate that, in principle, the accuracy of the numerical closure can be 
estimated by monitoring non-constrained state variables during the constrained simulation. 
Finally, when one can accurately assess the (lack of) accuracy of a given set of macroscopic 
state variables, it is straightforward to adjust the number of macroscopic state variables 
throughout a simulation using a corresponding accuracy criterion, as is done in [12, 18]. 
Note that, once a good set of macroscopic state variables is obtained, one could also consider 
proceeding along the lines of [19] to obtain a quasi-equilibrium closure. 

So far, we have not discussed potential gains in computational efficiency compared to 
a full microscopic simulation. One way to achieve a reduction in computational cost is 
to make use of coarse projective integration [24, 25] or similar methods [10, 11]. In this 
type of methods, one uses the proposed numerical closure technique to estimate the time 
derivative of the unavailable macroscopic model, and uses this estimated time derivative to 
perform a large (projective) forward Euler step for the macroscopic state variables; one then 
repeats the numerical closure procedure. The efficiency of coarse projective integration is 
strongly tied to a separation in time-scales between relaxation and macroscopic evolution; 
unfortunately, the physically interesting non-Newtonian behaviour precisely appears when 
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this time-scale separation is absent. We refer to [31] for a study of the acceleration that 
can be obtained in the small Deborah number limit, in which the polymeric fluid becomes 
Newtonian. 
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